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ABSTRACT 

A numerical model which employes observed stratospheric winds to 
advect simulated tracers was developed. This model successfully repro- 
duces many qualitative features of the observed fields of both ozone 
and radioactive debris. As the tracers evolve, the horizontal eddies 
constitute the principal process modifying the tracer zonal mean. 
Although the vertical eddies and the zonal mean cell are both an order 
of magnitude weaker than this process, the latter of these tends to 
always act in the same sense so that its effect becomes more important 
over long periods of time. 
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I . BACKGROUND 



A. OBSERVATIONS 

The study of trace substance transport in the stratosphere can 
provide many useful insights into the dynamics of that region. The 
history of an evolving inert tracer depends upon its initial state and 
upon several statistical properties of the wind field. In the strato- 
sphere, detailed observations of the wind field are not always plentiful 
enough to reveal these rather subtle characteristics of the general 
circulation. In addition to observations of the wind, there exist some 
observations of the distributions of such natural tracers as ozone and 
radioactive debris. These observations may be used to extend the avail- 
able data base, and theories of the general circulation may be tested 
by requiring them to simulate the advection of natural tracers. 

The influence of stratospheric transport processes was first noticed 
when the early spectrophotometer observations of Dobson and his co-workers 
(Dobson et al . 1928) revealed that the distribution of total atmospheric 
ozone has the following characteristics: 

1. The total amount of ozone per unit area increases from the equator 
toward the poles. 

2. In high latitudes this quantity varies with season exhibiting a 
pronounced maximum in late winter or early spring. 

3. No apparent correlation exists between total ozone and solar 
activity. 

4. Total ozone is correlated with the passage of upper-level synoptic 
systems, attaining positive departures from its long-term mean in the lows 
and negative departures in the highs. 
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Dobson's observations are consistent with those obtained by the 81 
station network in the Soviet Union and reported by Bozkov (1968). These 
observations indicate that: 

1. Minimum total ozone occurs at 10 N. 

2. Maximum total ozone occurs at 60 to 70 N, with values decreasing 
toward both the equator and the pole. 

3. The maximum zonal mean gradient of total ozone occurs between 30 
and 60 N. 

4. The total ozone tends to have a wavelike distribution in the 
horizontal with low values occurring over the western portions of conti- 
nents and high values over the eastern portions. 

The photochemical theory of ozone formation (See Craig, 1965) pre- 
dicts that, at photochemical equilibrium, the distribution of total 
ozone should attain a maximum in low latitudes at the subsolar point 
and increase toward the poles. Brewer (1949) and Dobson (1956) attempt 
to reconcile this theory with observation by postulating a mean cell with 
ascending motion in low latitudes, northward advection in the high mid- 
latitude stratosphere, descent near the pole, and a return circulation 
at low levels in mid-latitudes. This circularion supposedly carries the 
ozone, which is generated by solar ultraviolet radiation in low latitudes, 
northward and downward into the lower polar stratosphere where it ac- 
cumulates since there it is protected from photodissociation. 

The Brewer-Dobson theory is seriously challenged by airborne filter 
observation of the distribution of Tungsten-185 from the Hardtack nuclear 
tests (Feely and Spar, 1960). The center of concentration of this debris 
remained at 50 mb and nearly at the latitude of injection. The cloud did 
not migrate northward and downward, but rather spread in the meriodinal plane 
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in such a manner that the isopleths of concentration sloped down- 
ward from the equator toward the pole. 

B. PARAMETERIZED MODELS 

Reed (1950) and Normand (1953) note that stratospheric troughs are 
regions of both high temperature and high total ozone. This leads them 
to postulate that, since both potential temperature and ozone mixing 
ratio increase upward, the lows must be regions of descending motion 
and the highs regions of ascent. Furthermore, the troughs and ridges 
move more slowly than the wind so that horizontal advection shifts the 
centers of greatest mixing-ratio departure downwind from the associated 
pressure system. Because the wind almost always has a westerly component 
in the winter season, the ozone maxima lie in the southerly winds east 
of the troughs and the minima are east of the ridges in northerly winds. 
Thus, the departures of both temperature and mixing-ratio are explained 
in terms of the systematic interaction between horizontal and vertical 
advection. High mixing-ratio and temperature is accompanied by descent 
and poleward motion while low values of these two quantities occur with 
rising and equatorward motion. 

Citing earlier work by Dodson (1960), Newell (1961) carries this line 
of reasoning one step further. He explains the poleward slope of the 
isopleths by invoking eddy mixing in which northward motion is correlated 
with descent and southward motion with ascent. In the same paper he 
computes the flux of ozone by correlating total ozone with the winds at 
12-18 km in the maximum ozone layer. He finds that 90 percent of his 
computed flux is due to large-scale eddies and that the total flux is 
more than adequate to explain observed changes in the distribution. 
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Prabhakara (1961) generates zonal mean ozone distributions which 
depend upon the interaction of photochemistry, a mean cell that descends 
near the pole, and large-scale eddies parameterized as anisotropic dif- 
fusion. His results, which are in good qualitative agreement with obser- 
vation, show that such transport processes can explain many departures 
from photochemical equilibrium. 

In contrast to Prabhakara' s work, in which the principal axis of 
eddy diffusion is horizontal, there have been a number of models simu- 
lating the transport of radioactive debris with anisotropic diffusion 
schemes in which the principal axis is inclined to the horizontal (Reed 
and German, 1965; Davidson, et al . , 1965; Seitz, et al . , 1968; and Fair- 
hall and Reed, 1968). In such models it is assumed that the northward- 
downward velocity correlation in the large-scale eddies can be para- 
meterized in terms of a diffusion tensor whose principal axis slopes 
downward toward the pole. This means that the tracer is presented with 
a path of least resistance parallel to the observed slope of the iso- 
pi eths of mixing ratio in the real atmosphere. In order to keep the 
tracer on these sloping paths and to prevent excessive vertical spreading, 
it is necessary to make the vertical component of the diffusion tensor 

much smaller than that along the principal axis. Davidson et al . (1966) 

9 2 -1 

use a principal diffusivity of 4 x 10 cm sec and a vertical dif- 

4 2 -1 

fusivity of about 10 cm sec . However, since characteristic horizontal 

3 

distances are on the order of 10 km while vertical distances are about 
1 km and, since the speed of diffusion depends upon the square of these 

g 

distances, a ratio of 10 between the two diffusivities is not physically 
unreasonable. 
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C. GENERAL CIRCULATION MODELS 

Hunt and Manabe (1968) and Hunt (1969) report the most extensive and 
realistic simulation of stratospheric tracers to date. In their experi- 
ments various simulated substances are introduced into the wind fields 
generated by a general circulation model and the history of the numeri- 
cally evolved mixing ratio fields are examined. 

In the 1968 experiments inert tracers are permitted to evolve from 
zonally symmetric initial states. The first such tracer roughly corre- 
sponds to radioactive debris from an equatorial injection and is initially 
distributed in a band extending from the equator to ten degrees north with 
the maximum concentration at 50 mb pressure height. 

The second experiment's initial state represents the distribution of 
ozone at photochemical equilibrium throughout the region of integration. 
Both tracers are treated as inert substances inasmuch as, except for 
removal of any tracer that finds its way into the troposphere, no sources 
or sinks are included in these experiments. 

The simulated histories of both tracers are reassuringly similar to 
observations of real tracers in the atmosphere. Both the radioactive 
debris and the ozone begin at once to migrate northward and, to some 
extent, downward. Thus the concentration at the equator decreases and 
that at the poles increases, so that the initial gradient becomes smaller 
and eventually reverses. 

The actual transfer by the quasi-horizontal eddies is well illustrated 
by the evolution of the mixing ratio field at various horizontal levels in 
the model. The zonally symmetric character of the initial field quickly 
disappears becoming sinusoidal with wave number four predominating as it 
does in the circulation generated by the model. With the passage of time 
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the waves becomes increasingly exaggerated, forming the mixing-ratio 
contours into long northward protrusions. Eventually the protrusions 
separate, leaving isolated islands of high concentration in northern 
latitudes . 

Hunt (1969) reports two additional experiments with the same model 
and involving ozone transport with photochemical sources and sinks added. 
The mechanism of transfer in these studies and in those previously dis- 
cussed are very similar. In both, the mean cell, which produces a 
convergence of tracer into the subtropics, is opposed by an eddy diver- 
gence out of that region, with the eddies being particularly effective 
in transferring tracer northward in the region poleward of 30N. The 
principal difference introduced by photochemistry is the strong source 
in low latitudes which prevents reduction or reversal of the initial 
poleward gradient by replacing the ozone as fast as transport processes 
can remove it. 

The weak link in any experiment involving either parameterization 
or a general circulation model lies in the arbitrary specification of 
the transport process and its, at best, indirect relation to the actual 
motions in the stratosphere. In view of these difficulties, it seems 
instructive to compare advection of simulated tracers by actual observed 
winds with both observation and previous attempts to model both ozone 
and radioactive debris. The present study, therefore, employs actual 
observed stratospheric winds to advect simulated tracers in the presence 
of sub-grid-scale diffusion and numerically integrates the equation of 
continuity for the tracer in time to obtain the evolution of an inert 
tracer from an arbitrary initial state. 
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This approach is largely independent of any preconceived ideas about 
the stratospheric general circulation and so provides a third source of 
information intermediate between the parameterized models or the general 



circulation models and observations of real tracers. 



II. MODEL 



A. BASIC EQUATION 

The present study follows the lead of Hunt and Manabe (1968) by 
numerically exploring the detailed advection of trace substances with- 
out resort to arbitrarily defined diffusion tensors. Furthermore, the 
winds used to accomplish this advection are based on actual observations 
of the real atmosphere and are not the result of a mathematical model 
that may or may not reproduce the properties of the actual general 
circulation. 

The basic equation on which this study depends is the equation of 
continuity for an inert tracer without sources or sinks, expressed in 
flux form and spherical pressure coordinates. 



where m is the mixing ratio of the tracer, R the earth's radius, 0 the 
geographic latitude, /\ the longitude, p the pressure, CD the substantial 
pressure derivative following an individual air parcel, and u and v the 
wind's eastward and northward components respectively. 

The local mean of any dependent variable in (1) is defined as its 
average computed over a volume characteristic of the grid distance to be 
used in the integration. The value at any point within the region can 
now be represented by the sum of this mean and a departure from the mean 
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where the starred variables are the departures and the barred ones are 
the means. These relations are then substituted into (1) and the entire 
relation is averaged recognizing that, while the average of a fluctuation 
alone or of the product of a fluctuation with a mean vanish, the averaged 
product of two fluctuations does not. (See Haltiner & Martin, 1957). 



- od * i vcosy!) + u ^ 



The terms containing products of fluctuations are then approximated 
by using the Prandtl mixing-length assumption that the correlation between 
velocity and mixing-ratio is proportional to the mean gradient of mixing 
ratio as given by __ 

rfFcL* =■ -'K.Sffi 

K» 3rn (4) 



c)P 
no*u* - 



Rcos# 37 ^ 



The vertical and horizontal eddy diffusivities , and ^ are assumed 
to be constant in both space and time. This particular choice of con- 
stants, with the same eddy diffusivity in both horizontal directions and 
without any off-diagonal elements in the diffusivity tensor, implies that 
diffusion is isotropic in the horizontal plane with the axis of least 
diffusion oriented vertically since 

By substituting (4) for the eddy correlation terms in (3) and writing 
the flux terms back in advective form through use of the mass continuity 
equation for the local mean flow, one obtains 
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- - co 4- ^ QCO + u 
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— ^ 2|n0 _ — i-k 2. (r&bdi ^ ^£0 
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At this point (5) is placed in non-dimensional ized form by the fol 
lowing linear transformation of the velocity fields: 



(A =a 0 C\' j co - CS > 0 Co ' 



(6) 



Here the sub-zero values are constants with magnitudes (u Q = v Q = 20 kt, 
Ca) Q = 5 mb day - ') so chosen that the primed variables, whose variations 
contain the spatial and temporal changes of the velocities, are of order 
one. Since the equation is linear in m no such scaling needs to be done 
for thatvariable and (5) becomes 

- + VoV'om 0 \o uj ' dm 

c>t SP R_ Bp iz.cos0 

^ f —J ^5? *£25? (7) 

BPZ *(_ iZ^oSp D4>*- £Vos s <?$ 3VJ 



For solution on a digital computer the infinitesimal derivatives in 
(7) have to be replaced by finite differences, thus converting the partial 
differential equation into a system of algebraic equations in the mixing 
ratio and wind components at discrete points in the region of integration. 
This procedure requires that the independent variables of the problem be 
expressed as multiples of small but finite increments so that the co- 
ordinates of a point (p,G$,^, t) become i,j,k and n* where i = p/£ p, 
j = (P - > k = 1VSA and n = t/ &t and the values assumed by 

(p»0>?^ t) are required to make i, j, k and n integers. In this for- 
malism m(p,0,^, t) is ^ k and v(p, t) is v" 
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Once such a coordinate system is set up it becomes possible to 
approximate the time derivative by a forward difference 



gt > (8) 



the first order space derivatives by centered differences, for example, 
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and to use the three point approximation in place of the second order 
space derivatives 
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The /^operator always applies at the point i, j, k; although the space 
subscripts are usually omitted to simplify notation motion. The sub- 
script on the operator denotes the coordinate direction along which the 
difference is to be taken and the superscript the order of that dif- 



ference. 

With all the derivatives replaced by finite differences and after 
multiplication through by&t, (7) becomes 



rrT' = m° 

+ fS^l ^s^AziMAinrD f ~ K * 

(S )-l cost<4,- + uiVW’Jfo €&-&<})) 



(n) 



where the n and n + 1 superscripts denote the base time level and the 
newly generated time level respectively, and the £ superscript, which 
also denotes time, can take on only one or the other of these two values. 
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Each of the factors in brackets in (14) is the product of St 
with a collection of constants which, when combined, has the units of 
reciprocal time. Thus each bracket forms the dimensionless ratio of 
the timestep to a time characteristic of the particular process chang- 
ing the mixing ratio. For example, in the first term 2$p / q £ represents 
the vertical advection time or in the fourth term (%P) /K-j is the 
meridional diffusion time. 

Since each of these factors is a constant, (11) is simplified by 
replacing them with single quantities to reproduce the simpler form 
m n *' = m° - b,co J AliT)^- b a v^{m^-rb 3 /cosC<4-i 

~ K nrr*- Q bs /co <,(_ <jb 0 - £ <k><p)}\ (os(4> 0 - $ m* 

+ C bt /cos*(<A> - i A 2 *- fh ^ 

( 12 ) 



The pressure increment is 12.5 mb, that of latitude 5° and that of 
longitude 10°. With this information the first three advection times 
are evaluated to give: 

1\ - ^?/6Jo = b ciars ^ -U25 ddy^ 

V 3 = 2il^/u 0 = 2,5 days. ( 13 ) 

The value of does not represent the actual advection time enter- 
ing into the calculation since a factor of cos^S is not included in its 
evaluation. This factor represents the decrease in advection time as 
the distance corresponding to shrinks with approach to the pole. 
Since the cosine of the latitude varies from .819 at 35 N to .087 at 
85 N the actual zonal advection time changes by nearly an order of 
magnitude over the region of integration. 
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Since a state of nearly total ignorance exists with regard to the 
actual magnitudes of the eddy diffusivities in the stratosphere, the 
diffusion time along a given axis is set to a simple multiple of the 
advection time along the same axis. In this study the multiple (called 
1/b y) is chosen to be eight. This means that the parameterized dif- 
fusion can remove departures from equilibrium 1/8 as fast as advection 
operating by itself. This assumption leads to the following values for 
the diffusion times: 

^ - 4o d- ays , s 5 days 

1* - S>0 day5 . (14) 

The difference of a factor of four between Tj. and arises because 

= 2S^and, since the horizontal diffusion is assumed to be isotropic, 
the transport must occur four times as rapidly over half the distance. 

The values of and Kg corresponding to and are 1x10^ cm^ 
sec"^ and 4.5 x 10^ cm^ sec - ^ compared to 2 x 10^ and 1 x 10^ used in 
Prabhakara's parameterized study, That author's definition of diffusion 
and the one used here are essentially different. In the former case the 
diffusion was used to simulate transfer by the large-scale eddies, while 
in this study it has only to account for transfer by fluctuations too 
small to be represented in synoptic-scale data. The similarity between 
the two vertical diffusivities is unfortunate; one would have preferred 
to use a vertical diffusivity that differed from that of Prabhakara by 
a factor of at least ten because, in the parameterized studies, diffusion 
is called upon to provide much of the transport that advection accomplishes 
in this model . 

The entire treatment of diffusion is rather artificial but since it 
serves as a computational smoother, its effect can not be minimized with- 
out sacrificing numerical stability. 
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B. TIME INTEGRATION 



Integration in time is accomplished using a backward corrected 
Euler scheme (see Kurihara, 1965). If the advective and diffusive 
terms in (12) may be replaced by f(m n , n , v n , u n ) this scheme may be 
represented thus: 

m £ = + f OT)° } co 1 ^ v'] CAP ) 



w 0H = m°4. (cmi ^ ; v 0+ ; u n ") m 



(15) 

(16) 



Given the field of m n and the wind at the n time level tn^ is computed 
from (15), Then m A and the winds at the n + 1 time level are used in 
(16) to obtain m n+ b Thus m n ^ is first predicted using a forward 
difference and then corrected with a backward difference, the entire 
process being repeated at each successive time level until m is known 
for the entire period of interest. 

According to Kurihara (1965), if this method is applied to the one 
dimensional advection equation, 

- 9(°-h b A; X' 1 



-X mi r bA i *** 



07) 

(18) 



stability will result for b ^ 1 . By experiment, it was discovered that 
optimum results are obtained by using a time step £>t = 1/48 day producing 
the following values for the b’s (b^ = 1/8) 

b, -4.i7*id 3 , /. 67 >* 10 ® 63 -^. 33^10 3 

b* = o. 5 ?x 10 b 5 - 4 . 17* 10 i j b t -).o 4 *io 3 (19) 

These values are obviously much less than one, but the factor of 1/cos^ 
makes the coefficient of the third term much larger near the pole. At 
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80 N, the most northerly point at which (15) and (16) are applied, 
cos = .174 so b^/cos - 0.047 still well within the region of 
stability. Unfortunately , a longer time-step drastically reduces the 
quality of the results and, although no theoretical investigation is 
attempted, this may be caused by either mass imbalances in the wind 
data or by the variability of the space increment resulting from the 
cos^ effect. 

C. WIND DATA 

The wind data used by Mahlman (1967, 1969) are employed in the 
advective terms of (11). The data, which cover the 41 day period from 
15 November until 25 December 1958, consist of the three wind components 
and the fields of temperature and geopotential height tabulated at 50 
and 100 mb for the region from 40 to 80 N. Originally the two hori- 
zontal wind components as well as the temperature and height fields 
were extracted directly from the Weather Bureau (1963) stratospheric 
daily chart series for the IGY, while theca's were computed from the 
other fields using the thermodynamic equation with an assumed diabatic 
heating rate (Mahlman, 1967). 

During the early portion of the period of data coverage, the polar 
night vortex is intensifying with the circulation dominated by a cold 
low over northern Asia. Then on 21 November the flow undergoes a pulsa- 
tion or "minor breakdown", becoming more meridional in character and 
remaining disturbed until 11 December when the vortex re-stabil izes with 
an important wave number two component. During this last period one 
trough lies over Asia and the other over North America with the two 
intervening ridges over the oceans. 
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With the exception of the ten-day period from 5 December until 15 
December, the mean cell in these data has rising motion over the pole 
and descent in middle latitudes (Mahlman, 1969). Even during the time 
in December when the vertical motion reverses at the pole, descent takes 
place only in the immediate region of the pole with rising motion 
continuing in the 60 to 70 N band. This means that, except for that 
brief period in December, the mean cell oeprates in the exact opposite 
sense to that required by the Brewer-Dobson theory. 

Unfortunately , the data are tabulated every 24 hours, on only two 
pressure surfaces, and for a somewhat coarser horizontal grid. This 
means that the data have to be extended by extrapolation and interpolation 
in both space and time- First, using the thermal wind equation, the hor- 
izontal winds are extrapolated upward to 12.5 mb and downward to 150 mb, 
and the <L>'s at these two levels are assumed to be zero. This is done 
at every point in the horizontal that carries the initial data fields 
(every 40° of longitude at 80 N, every 20° at 70 N, and every 10° of 
longitude at 50, 50 and 40 N). Values for every 10° longitude at 70 N 
80 N are generated by linear interpolation between the tabulated values 
and the same process is applied again to produce winds at odd multiples 
of five degrees of latitude. Finally, a cubic Lagrange interpolating 
polynomial is fitted to the data at each horizontal grid intersection for 
the four levels and evaluated at the intermediate pressure levels to 
complete the wind field for each day The wind for each time step is 
obtained from the daily wind fields by simple linear interpolation in 
time. 

The wind data require still more treatment before they are suitable 
for use in the integration scheme- Due to inaccuracies in observation 
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and processing, the probable errors in the v wind component are at least 
as large as the real mean cell. Furthermore, the omegas generated by 
the thermodynamic method are not necessarily consistent with the diver- 
gences of the horizontal wind. Mass imbalances in the large-scale eddies, 
while undesirable, tend to cancel each other when integrated over the 
entire region, but a mean cell with spurious mass sources may seriously 
distort the results of the computation. It was decided to remedy the 
situation by removing the natural mean cell and replacing it with an 
artificial time independent one which exactly satisfied the zonally- 
averaged continuity equation. The following transformation accomplishes 
this end: 




V 



i) A 



h 







- CO 



"A 

‘‘is 



- 



+n-t 



s 



( 20 ) 



Here the double primed quantities are the new values, the unprimed the 
old, the tilded the zonal mean old values, and the capitals the contrived 
mean. This new mean is usually run with ascent at the poles; although 
some experiments were run with the cell reversed for comparison. Typical 
magnitudes for the new mean were on the order of a knot for V and .5 mb 
day - ^ forX2.. The new cell was constant in time and no attempt was made 
to introduce time variations resembling those of the natural mean cell. 

While use of an artificial mean cell seems to introduce the very element 
of arbitrariness that was avoided through employing observed winds, it 
should be pointed out that this approach both permits experimentation with 
mean cells of differing character and allows use of mean circulations that 
are more or less consistent with heat and momentum considerations. 
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Ill, EXPERIMENTS 



A. OZONE SIMULATION 

1 . Boundary and Initial Conditions 

The behavior of stratospheric ozone is simulated by the evolution 
of an inert tracer from a zonally symmetric initial state roughly cor- 
responding to the observed zonal mean ozone distribution (Hering, 1960). 
Throughout the integration the observed initial values at 12.5 and 150 
mb are retained at the top and bottom boundary points. This condition 
can be justified since it is assumed that the value at 12.5 mb is control- 
led only by photochemistry and is time-independent except for diurnal and 
seasonal variation which are not reproduced. Similarly the mixing-ratio 
at 150 mb is assumed to represent a constant equilibrium between down- 
ward transport and destruction in the troposphere. 

Originally it was the intention to parameterize advection into 
the region from the tropics by retaining the observed zonal mean values 
at the southern boundary points, but the final version of the model is 
permitted to generate its own southern boundary condition inasmuch as 
the points at 35 N are all set to the generated zonal mean of those at 
40 N. 

Near the pole, the cosine of the latitude becomes small causing 
an apparent singularity in the prognostic equation. This difficulty can 
be avoided by applying a cartesian version of the prognostic equation at 
the pole, but a simplier solution is suggested by the treatment of the 
southern boundary, so the values at 85 N are set to the zonal mean of 
those at 80 N, 
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These boundary conditions imply that the system is not closed 
since boundary fluxes can occur through all surfaces of the region of 
integration. A horizontal flux through the southern boundary is physically 
reasonable since advection from tropical regions is an important source 
of mid and high-latitude ozone. Similarly, the flux over the pole is 
to be expected; although the boundary scheme used does not properly 
simulate it, and a spurious over-the-pole transfer arises which distorts 
the final results. The upper boundary condition allows ozone to be 
transported downward into the region of integration but does not permit 
it to leave through the top. The values at 12,5 mb are invariably higher 
than those in the interior of the region. Since the vertical motion is 
small near the upper boundary, diffusion predominates allowing the tracer 
to flow inward but never outward. A similar set of conditions arises at 
the lower boundary; except here the boundary values are always lower than 
those in the interior so only outward fluxes can occur. Thus, although 
no source of sink terms are explicitly included in the model, the boundary 
conditions crudely approximate the behavior of photochemical sources in 
the upper stratosphere with a sink at the tropopause. 

2. Results 

Figures 1 through 5 are a series of charts portraying the ozone 
mixing ratio and geopotential height on the 50 and 100 mb levels at 
intervals of 5 days during the period of integration. Comparison between 
the 50 mb and 100 mb pressure levels indicate that the degree of continuity 
of the ozone field between these two levels is comparable to that observed 
in the geopotential height field so that representation of more than one 
horizontal section of the ozone field is superfluous, and only the 100 mb 
surfaces are presented after 25 November, 
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lines are jeopotential height contours. 
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FIG. 2. Simulated ozone mixing ratio on 25 November at A: 





31 



FIG. 3. Simulated ozone mixing ratio at 100 mb on A: 30 November; and B: 5 December. 
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FIG. 4. Simulated ozone mixing ratio at 100 mb on A: 10 December; and B: 15 December. 
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FIG. 5. Simulated ozone mixing ratio at 100 mb on A: 20 December; and B: 25 December. 
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FIG. 6. Simulated ozone mixing ratio on 20 December for A: reverse mean cell; 

and B: zero mean cell. 






December with mean ascent at the pole; C: for the 

same date with a zero mean cell; and D: for a mean 

cell with descent at the pole. 
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These charts show a number of qualitative features that are 
reassuringly similar to the characteristics of the real atmospheric 
ozone field. Dobson's (1924) observations of ozone maxima in the 
troughs are reproduced. Similarly, Bozkov's (1968) observations of 
sinusoidal isopleths of concentration with ozone maxima in the standing 
troughs over continents and minima over the intervening seas appear in 
these fields. 

At 50 mb and 60 N the normalized correlation coefficient between 
o and m, computed over the entire period of integration, is 0.22. This 
means that the high mixing ratios in the lows can be explained in terms 
of advection from above just as Reed (1950) and Normand (1953) postulated. 

On the 50 mb surface the ozone mixing ratio values initially 
decrease from .053/igm/gm at 35 N to 0.028 / 4<gm/gm at 85 N while at 100 
mb the initial values increase from 0.010 to .031 in the same interval. 

As the integration proceeds the transport processes reverse the 50 mb 
gradient by 18 November and establish a mid-latitude maximum at 75 N 
by 25 November. This maximum tends to drift southward appearing most 
frequently at about 65 N during the latter portion of the experiment. It 
is quite persistent in this latitude, but often it disappears for several 
days, being absent on a total of 19 of the 41 simulated days. A similar 
maximum appears in the 100 mb profile on 2 December but it is less pro- 
nounced and not nearly so persistent in time (being present on 10 of 41 
days). At both altitudes the maximum zonal mean gradient occurs south 
of the latitude of maximum mixing ratio, with the field being more nearly 
flat over the pole. These results are consistant with Bozkov's work, and 
since they arise during the integration and do not represent more persist- 
ence of the initial state, they enhance confidence in the model. 
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When time series of the zonal mean mixing ratio and of the 
various processes modifying it are considered (see FIG. 8), it is found 
that the horizontal eddies and horizontal diffusion are very highly cor- 
related with the fluctuations of the mixing ratio. The other processes, 
the vertical eddies, the vertical and horizontal mean cells and the 
vertical diffusion are all about an order of magnitude smaller than 
these first two effects and did not produce noti cable, short-term changes 
in the mixing ratio. 

The direction of horizontal eddy transport changes during the 
integration. At the beginning it is poleward, but during the pulsation 
of the polar night vortex it changes sign and remains negative for 22 
days until the vortex begins to re-intensify. When the normalized mixing- 
ratio and meridional velocity correlation coefficient is computed for the 
entire period, its value is + .023 at 50 mb and 60 N. But when the data 
are stratified into poleward and equatorward transport periods, the cor- 
relation during the former is + .317 and the latter -.237. Since the 
omega-mixing-ratio correlation is + 0.22, this indicates that the dif- 
ference in the effectiveness between the vertical and horizontal eddies 
is primarily caused by the smaller ratio of characteristic velocity to 
space increment for the former. 

Because the mean cell is constant in time and the zonal mean 
gradient is nearly so, the sign and magnitude of the instantaneous change 
resulting from the mean cell is quite persistent. This means that, while 
the short-term effect is small, the mean cell produces a marked change in 
the field when integrated over 41 days. 

This effect is well illustrated by the differences among the 20 
December 100 mb fields in FIG. 5 and FIG. 6. Figure 6 shows the 20 



38 



December mixing-ratio field resulting from a zero mean cell as well 
as that resulting from a reversed mean cell. Comparison among these 
three charts indicates that in all cases the features are very similar 
in shape and location but the values differ as though the fields were 
translated by addition of a constant. This is very nearly what happened. 
Figure 7 shows the zonal mean cross sections for each of these three 
cases on 25 December as well as the initial zonal mean. Each different 
mean cell changes the inclination of the mixing ratio isopleths, tending 
to increase their slopes when sinking occurs at the pole and decrease 
them when rising takes place at the pole. This effect is most important 
in the northern portions of the region where it moves the lines in the 
vertical while toward south they remain at about the same altitude in 
all cases. It should be pointed out that, while this effect becomes 
important only after long integration times, it will eventually attain 
an equilibrium with the other processes and not accumulate indefinitely. 

Apparently the model seeks to attain just such an equilibrium 
during this simulation. At 50 mb and 60 N the mixing ratio rises 
asymptotically to a value of about 60 yUgm/gm”^ (See FIG. 8). The 
horizontal eddies, the vertical diffusion, and the horizontal mean all 
appear to contain decaying transients that largely damp out by 10 December. 
Thus, apparently about 25 days are required to recover from the dis- 
equilibrium introduced when the model is initiated from a zonally sym- 
metric initial state. 

The contribution of both components of the sub-grid-scale diffusion 
is disturbingly large- This phenomenon is dependent upon a completely 
arbitrary sub-grid-scale parameterization scheme and probably is much 
stronger than any real atmospheric diffusion. It is unfortunate that 
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the diffusivities can not be reduced without impairing the model's 
computational properties. The un-naturally strong diffusion may mask 
important advective processes and spuriously shorten the equilibrium 
time as well as displace the equilibrium from its true value. 

B. RADIOACTIVE DEBRIS SIMULATION 

1 . Boundary and Initial Conditions 

Clouds of radioactive debris are simulated by permitting the 
distribution of tracer to evolve from an initial three dimensional 
Gaussian distribution centered at a point within the region of integra- 
tion. It would be more satisfying to use an initial state in which all 
the tracer is concentrated at a single grid point, but in such experi- 
ments an undesirable computational mode is excited by the steep gradients 
so that the results are meaningless. 

The boundary conditions imposed upon the evolving cloud of debris 
represent a serious problem. For the previous experiment with zonally 
symmetric ozone fields as an initial state, the generated zonal mean of 
points along the next interior latitude circle represents a reasonable 
boundary value. Because the point source problem is asymmetrical, con- 
straining the gradient at the southern boundary to zero seems to constitute 
a reasonable boundary condition, but computational difficulties arise in 
this case when the wind normal to the boundary is strong. The final 
solution to this problem lies in a compromise between the two approaches. 
Each boundary point at 35 N or 85 N is set to a weighted average of the 
zonal mean at the next interior parallel and the value at the adjacent 
point along that parallel. The weighting factors are 0.7 times the mean 
and 0.3 times the adjacent point at the southern boundary^ and equal weights 
are used at the northern boundary. 



40 




41 



FIG. 9. Simulated radioactive debris distribution at 100 mb in arbitrary units on 
18 November; and B: 20 November. The heavy lines are the mixing ratio 

contours while the light lines are the geopotential height contours. 
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FIG. 10. Simulated radioactive debris at 100 mb on A: 22 November; and B: 24 November. 








FIG. 11. Zonal-mean simulated radioactive debris in arbitrary 
units on A: 18 November; B: 20 November; C: 22 

November; D: 24 November; and D: 26 November. 
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At the top and bottom of the region there is little inflow or 
outflow because the vertical velocities are small. Furthermore, these 
regions are fairly remote from the centers of concentration. This means 
that zero gradient at 12.5 and at 150 mb constitutes a physically reason- 
able and numerically stable boundary condition there. 

2. Results 

The evolution of a cloud injected at latitude 60 N longitude 
140 W with initial standard deviation 10° of latitude, 20° of longitude 
and 25 mb of pressure was simulated for eleven days (See FIGS. 9, 10). 

This initial position lies in the almost purely zonal flow in the ridge 
over the Northern Pacific so that the main mass of debris is carried 
eastward along the height contours while gradually spreading and dis- 
persing. By 23 November the cloud has migrated about 90° of longitude 
and lies in the eastern portion of the Atlantic ridge. During this period 
the mixing ratio at the last closed contour decreases from 80, its initial 
value, to 20 (arbitrary units). 

On 17 November part of the cloud is carried over the pole and 
then southeastward around the Siberian low. Between 19 and 21 November 
the mass rounds the trough and begins to move northward so that by 23 
November it has circumnavigated the pole at high latitudes and re-merged 
with the main cloud. 

In the zonal mean, the center of the initial distribution lay at 
60 N and a pressure altitude of 75 mb (See FIG. 11). As the cloud evolves 
the center of maximum mixing ratio gains altitude and moves toward the 
south, at the same time spreading laterally in such a manner that its 
axis slopes downward toward the pole. By 23 November the zonal mean 
bears a striking resemblance to both the results of parameterized studies 
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and to observations (Gudiksen et al . , 1969), showing that this model can, 
through advection and diffusion alone, quickly shape any reasonable dis- 
tribution of inert tracer into a fair reproduction of observed mixing 
ratio fields. In fact, this process took place much more quickly than 
it would have in the real atmosphere because of the large standard 
deviation in the initial distribution. 

No detailed analysis is made of the principal agencies responsible 
for the tracer's evolution, but the horizontal eddies and diffusion com- 
bined are about an order of magnitude stronger than the sum of the other 
four agencies. However, as in the ozone experiments, the sign of the 
tendencies produced by the mean cell terms remain consistent in time so 
that over periods longer than eleven days their effect becomes significant. 

IV. CONCLUSIONS 



This model successfully simulates many observed qualitative features 
of stratospheric transports of both ozone and radioactive debris. Apart 
from some numerical difficulties which result from mass imbalance in the 
wind field and from inadequate treatment of the singularity at the pole, 
no obstacles exist to prevent the use of this model to predict the evolu- 
tion of stratospheric tracers. 

If the erroneous exaggeration of the parameterized diffusion may be 
neglected, the instantaneous derivative is primarily determined by the 
horizontal eddies. Although the effect of the vertical eddies is about 
an order of magnitude weaker than that of the horizontal eddies, it can- 
not be neglected in the zonal mean. The two components of the mean cell 
are still weaker than the vertical eddies. However, over long periods of 
time, they may assume greater importance because they do not change sign 
or fluctuate greatly in magnitude. 
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